Glass Transition of Hard Sphere Systems: Molecular Dynamics and Density 

Functional Theory 



m 
O 
O 

(N 

Oh! 
< 
(N 



Kang Kim and Toyonori Munakata 
Department of Applied Mathematics and Physics, Graduate School of Informatics, 
Kyoto University, Kyoto 606-8501, Japan 
(February 1, 2008) 

The glass transition of a hard sphere system is investigated within the framework of the density 
functional theory (DFT). Molecular dynamics (MD) simulations are performed to study dynamical 
behavior of the system on the one hand and to provide the data to produce the density field for the 
DFT on the other hand. Energy landscape analysis based on the DFT shows that there appears a 
metastable (local) free energy minimum representing an amorphous state as the density is increased. 
This state turns out to become stable, compared with the uniform liquid, at some density, around 
which we also observe sharp slowing down of the a relaxation in MD simulations. 
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Understanding of the universal mechanism of the glass 
transition is one of the major challenges in the current 
condensed matter physics. From a dynamical viewpoint, 
we would like to know how a drastic slowing down near 
the transition point (temperature or density) and an 
eventual exceeding of the relaxation time over the ex- 
perimental time scale are realized [1,2]. Energetically or 
statically, it is asked whether a thermodynamic glassy 
state with the amorphous arrangement of particles has 
lower free energy than a liquid state of uniform density. 
In order to answer these questions, many efforts have 
been devoted to real experiments, computer simulations, 
and theories for the last decades. 

As one of theories to study supercooled liquids and 
glasses, the density functional theory (DFT) is recently 
gathering much attention [3]. The DFT is now a con- 
ventional method to study the freezing [4,5] and other 
transitions. The glass transition has been investigated 
also based on the DFT by some workers [6-10]. In the 
earlier works [6-10], the random close packing (RCP) of 
hard spheres has been produced by the Bennett's algo- 
rithm [11] and the free energy from the DFT with the 
input density field supplied by the RCP data has been 
calculated. Singh et al. [6] showed that the hard sphere 
glassy state becomes more stable than a uniform liquid at 
the critical density n g a 3 = 1.14, suggesting that there ex- 
ists a kind of the thermodynamic (later called a random 
first order [12]) glass transition. Here a is the hard sphere 
diameter and n is the number density of the system. It 
is remarked here that since the RCP configurations were 
produced by a kind of aggregation method, we can not 
study dynamical aspects of the glass transition found by 
the energetics based on DFT. 

The purpose of this paper is first to produce a su- 
percooled and a glassy state for a one-component hard 
sphere system, relying on the uniform compression molec- 
ular dynamics (MD) method recently developed by 
Lubachevsky and Stillinger [13] and then to study both a 
dynamical and static properties of the state. Especially 



from the particle configuration data, we can discuss free 
energy within the DFT framework. This MD approach, 
in conjunction with the DFT, enables us to study both 
dynamical and static aspects of the transition in contrast 
with the Bennett's approach. 

Our system consists of N = 1372 identical hard spheres 
with the mass m and the diameter a in a cubic box of 
volume V with periodic boundary conditions. Through- 
out this paper, the units of length and time are a and 
yjmcr 2 /UbT , respectively, where fee is Boltzmann's con- 
stant and T is the temperature [14]. It should be noted 
that the temperature is fixed as kgT = 1 in the course 
of MD simulations. 

To begin with, we briefly explain our MD method to 
obtain glassy states of a one-component hard sphere sys- 
tem. Employing the standard Alder and Wainwright 
algorithm [15,16], we first generate the equilibrium liq- 
uid state at the density n = 0.86. It is well known 
that the fluid system freezes at n/ ~ 0.94. To avoid 
the crystallization and obtain amorphous glassy states, 
Lubachevsky and Stillinger introduced a compressing (or 
quenching) procedure [13], in which they actually in- 
creased the diameter a with a constant rate of expansion 
during MD simulations. The dimensionless expanding 
rate T is defined as 
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and r = 0.01 is chosen in our simulations. From the ini- 
tial state h — 0.86, we expanded each sphere with the 
rate T and could obtain various high density states, 0.86, 
0.94, 1.02, 1.06, 1.10, 1.14, 1.18 and 1.21, without crys- 
tallization. 

Let us first study the static structure of the system. 
For the purpose the radial distribution function g(r), 
which is denned by 
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is calculated, where v% represents the positions of the i- 
th particle and (• • • } denotes the ensemble average over 
different configurations. In Fig. 1, we plotted g(r) for 
the density n = 0.94, 1.06, 1.14, and 1.21. It is noted 
that there is no sign of crystallization, which would be 
reflected in sharp peaks of g(r) at some characteristic lat- 
tice spacings. Instead, g(r) for higher density (h > 1.14) 
shows the splitting of the second peak, which is a fa- 
miliar characteristic of a glassy state of a simple liquid. 
Furthermore, we notice that the contact value g(r = a) 
shows an anomalous increase as the density increases (see 
the inset), which corresponds to the fact that the pres- 
sure increases drastically with increasing n. Incidentally 
it should be remarked that the forms of g(r) agree qual- 
itatively with those illustrated in Rcf. [11]. 

We next consider the dynamical aspects of the hard 
sphere glasses by calculating the incoherent intermediate 
scattering function F s (q,t), which is defined by 
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where Arj(t;to) = rj(t + t ) — rj(t ) is the displace- 
ment vector of the j-th particle in time t, and (• • • }t 
presents an average over initial times to- It is noted that 
the F s (q,t) is one of the standard quantities in stud- 
ies of dynamical properties of supercooled liquids and 
glasses [17,18]. In Fig. 2, we plotted the decay profiles 
of F s (q, t) at a dimensionless wave number q = 2tt for 
n = 0.94, f .06, 1.14, and 1.21. We see in Fig. 2 that the 
relaxation of the F s (q,t) at n = 0.94 can be expressed 
by a simple exponential function. Beyond the density 
h = 1.06, however, F s (q,t) exhibits the two-step, that is 
fast (3 and slow a, relaxation, which is often mentioned as 
a characteristic sign of the slow relaxation in glass form- 
ing liquids. At the highest density n = 1.21, the F s (q,t) 
dose not show any decaying behavior [19]. One can define 
the structural relaxation time f by F s (q = 2tt,t) = er 1 
and this f is plotted as a function of n in Fig. 3. We 
notice in Fig. 3 that the relaxation time f shows a strong 
dependence on the density, which can be expressed by 
the power-law (n Sl MD — n) -7 with n g ^MD — 1-15 and 
7 ~ 1.31 (solid line in Fig. 3). 

We now consider energetics of the system based on the 
DFT and the configuration generated by the MD simula- 
tions. For a practical calculation of the DFT, we employ 
the Ramakrishnan and Yussouff (RY) free energy func- 
tional [4] because of its simplicity and physical clarity. 
The RY functional F[n] is given by 



F[n]=F id + Fi$+F, 



where 
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F^ = ~k B Tj J[n(r)-n]C(\r-r'\) 

x [n(r') - n]drdr'. (6) 

Here, F^ and F uni represent the ideal gas contribution 
and the excess free energy of the uniform liquid state 

( 2) 

n(r) = n, respectively. F> n j. represents the second order 
term in the expansion around the uniform liquid state, 
thus all terms higher than second neglected. We note 
that C(r) is the direct correlation function of the uni- 
form liquid with the density n [20]. 

In order to evaluate the free energy of the system, we 
need the trial density field n(r), for which we employ a 
conventional Gaussian superposition [6-10]. That is, the 
density field n(r) is expressed by a sum of Gaussians with 
the centers located at N sites {ri}, which are given by 
our MD simulation. 

N ,,3/2 N 

n ( r ) = ^2[-) exp[-a(r -r 4 2 )] = ^z(r;ri), (7) 

i— 1 ^ ^ i— 1 

where a -1 , a variational parameter for the calculation of 
the free energy, is proportional to the mean square dis- 
placement of each particle. Small (large) a represents 
the uniform liquid (localized amorphous) state. 

When a is very large, F^ is asymptotically represented 

by [6] 
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For a small a region, we calculated the integral Eq. (5) 
numerically. We confirmed that F id approaches when 
<5^0 and noticed that Fid coincides with the value of 
Eq. (8) for a > 20. 

It is easy to see that the interaction term Fi nt can be 
divided into three parts as [6] 



^1 2 *» = ~Nk B T{ F inttl (a) + F intt2 (a) 
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where F int .i(a) represents the self-interaction of a single 
Gaussian, 

F int ,i(a)=J J z(r;0)C(\r-r'\)z(r':0)drdr' 1 (10) 

and -Fmt, 2 (a) represents the interaction between the two 
distinct Gaussians, 

F int , 2 (a)=n Jg(n)dn J J z(r;0)C(\r - r'\) 

xz(r';ri)drdr' . (11) 

The pair distribution function g(r) in this equation is 
determined from the MD simulation. With respect to 
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the direct correlation function C(r), we use Henderson- 
Grundke expression for C(r), which is known to be reli- 
able, though empirical, even for high density hard sphere 
liquid [21]. 

The total free energy per particle relative to uniform 
state, 



A/(a) = 



F id (a) + F%(a) 
Nk B T 



(12) 



is calculated as a function of the localization parameter 
a (see Figs. 4). It is seen in Figs. 4 that the free energy 
local minimum at finite a, which represents an amor- 
phous state, appears as the density is increased. As is 
mentioned in Ref. [6], the local minimum appears as the 
result of the competition between the ideal gas (simple 
increasing function of a) and the interaction terms. Fig- 
ures 4 show that two local minima are located at a ~ 13 
and 1600 for h = 1.14. Das et al. also observe two 
local minima of A/(a), which are called the weakly lo- 
calized state for small a and the highly localized state for 
large a, respectively [9]. In addition, similar values for a 
are reported in Refs. [6,9], which also use the RY func- 
tional, in relation to the local minimum of A/ (a). As 
is stated in Ref. [9], the qualitative adequacy is ambigu- 
ous for the highly localized state with very high value 
of a since the RY form includes a perturbative expan- 
sion around the uniform state. In fact, based on the 
MD data we estimated a for high density states by relat- 
ing it to the plateau value of the time dependent mean 
square displacement of each particle, yielding a ~ 50 for 
n = 1.14 for instance. From this it is seen that the RY 
form does not give the proper estimation of the degree of 
localization compared with the results obtained in earlier 
works [8,10]. 

In Fig. 5, we plotted the free energy differences A/ 
of the weakly and highly localized states as a function 
of the density h. From Fig. 5, we find that both the 
weakly and highly localized states become more stable 
than the uniform state at n g .DFT — 1.15, which is the 
liquid-glass transition density from the energetics based 
on the DFT. Moreover, it is seen that the highly localized 
state is more stable than the weakly localized state for 
higher densities. In passing we note that our (first order) 
glass transition density hg t BFT = 1-15 is rather close to 
the one h g = 1.14 found in Ref. [6]. 

Finally, we compare our results from the energetics 
above with dynamical informations supplied by our MD. 
We find in Fig. 2 that the intermediate scattering func- 
tion F s (q, t) begins to exhibit the two-step relaxation at 
the density n ~ 1.06, which corresponds precisely to the 
density where the free energy local minimum begins to 
appear in our DFT (see Figs. 4). Turning to the relax- 
ation time f , we recall that the density dependence of f 
could be described by the power-law {n gt MD — n)^ 1 with 
n g> MD — 1-15. This density happened to coincide with 
the density n g ,D ft beyond which the localized state is 
more stable than the uniform liquid in the present DFT. 



From these results, we consider that the DFT based ener- 
getics and dynamical behaviors related to slow dynamics 
are well correlated with each other. 

In this paper, we reconsidered the DFT approach to 
the glass transition in the hard sphere system, which was 
first undertaken by Singh et al. We obtained hard sphere 
glasses by MD simulations without recourse to the Ben- 
nett algorithm and the information on particle config- 
urations produced by the MD simulations are used as 
input data when the free energy is calculated based on 
the DFT. While only the uniform liquid state is stable at 
low density, the free energy local minimum begins to ap- 
pear at high density n ~ 1.06, where our MD shows that 
two-step relaxation begins to appear. This metastable 
glassy state becomes stable relative to uniform liquid at 
n g ,MD = 1.15. Slow relaxation as represented by the 
F s (q,t) turned out to be consistent with the energetics 
based on the DFT. 

Before concluding this paper, we comment on recent 
developments in studies of the DFT. In recent years, the 
so-called weighted density approximation (WDA) for the 
free energy functional has been developed [22,23] and the 
modified version is also introduced [24]. Moreover, the 
fundamental measure theory (FMT) is proposed [25] and 
gathering considerable interest. It is well known that for 
highly localized states, such methods are more accurate 
than the RY one. This is because the former employs a 
non-perturbation approximation whereas the latter relies 
a perturbation expansion around the uniform state. Sev- 
eral workers have already investigated the glass transition 
by using the modified WDA method and found that the 
metastable localized state is located at a ~ 100 [8,10] 
in accordance also with our MD results . Although the 
DFT based on the RY functional is still useful because 
of its physical clarity, generality and simplicity, we think 
that in view of the recent achievements it is meaningful 
to employ the WDA or FMT method in order to obtain 
improved results. 

Furthermore, we hope that the present DFT approach 
will be applied to more complex systems. As a model of 
a glass forming liquid, binary Lennard- Jones [17] or soft- 
core system [18] has been investigated by large scale MD 
simulations. Our approach can be readily applied to such 
system and would give new insights into the glass tran- 
sition from both thermodynamic and dynamical view- 
points. 
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FIG. 1. The radial distribution function g(r) obtained for 
h — 0.94 (solid line), 1.06 (dashed line), 1.14 (short dashed 
line), and 1.21 (dot dashed line). Inset: contact value of radial 
distribution function g(r = a) as a function of n. The units 
of f and h are a (f = ra) and <r~ 3 (n = na 3 ), respectively. 
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FIG. 2. Intermediate scattering function F s (q,t) at a wave 
number q = 2ir for n = 0.94 (solid line), 1.06 (dashed line), 
1.14 (short dashed line), and 1.21 (dot dashed line). The units 
of q and t are a (q = qa) and \J ma 2 jksT (i = t\JksT /ma 2 ), 
respectively. 
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FIG. 3. Structural relaxation time f as a function of den- 
sity n (closed circles). Solid line represents power-law fit 
(n g ,MD — n)~ 7 with h g .MD = 1.15 and 7 = 1.31. 



FIG. 4. Total free energy per particle relative to uniform 
liquid A/(a) for a < 100 (a) and a > 100 (b) as a function 
of localization parameter a. The unit of a is ct~ 2 (q = aa 2 ). 
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FIG. 5. Free energy differences A/ of the weakly localize 
state (solid line) and the highly localized state(dashed line) 
as a function of density n. 
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